Nonequilibrium Green's function method for thermal transport in junctions 
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We present a detailed treatment of the nonequilibrium Green's function method for thermal 
transport due to atomic vibrations in nanostructures. Some of the key equations, such as self-energy 
and conductance with nonlinear effect, are derived. A self-consistent mean-field theory is proposed. 
Computational procedures are discussed. The method is applied to a number of systems including 
one-dimensional chains, a benzene ring junction, and carbon nanotubes. Mean-field calculations of 
the Fermi-Pasta-Ulam model are compared with classical molecular dynamics simulations. We find 
that nonlinearity suppresses thermal transport even at moderately high temperatures. 
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I. INTRODUCTION 

Fourier's law describes the transport of heat in a 
macroscopic object. The calculation of the thermal con- 
ductivity that appears in the Fourier's law is fundamental 
and important in understanding the properties of materi- 
als. Such a calculation must be grounded upon quantum 
mechanics with phonons as basic quasi-particles. This is 
a non-trivial task and was successfully carried out many 
years ago by Peierls [3, 0] using the Boltzmann equation 
for phonons. Molecular dynamic (MD) simulation is an- 
other approach. However, MD results are not correct at 
low temperatures as it is purely a classical approach. 

The early works are mostly concerned with bulk ma- 
terials with many degrees of freedom and periodic lat- 
tices |3|. In recent years, more attentions are paid to 
heat transport in small or low-dimensional systems such 
as carbon nanotubes and molecular junctions 0, H[. In 
nano- to meso-scales, new features come in. In most of 
these conditions, the discreteness of the atom is impor- 
tant. The concept of phonons developed for the periodic 
lattices is somewhat difficult to apply if a system does not 
possess translational invariance, such as a nano-junction. 
In such situation, the Boltzmann-Peierls' approach is not 
applicable. 

There are several alternative methods that have over- 
come the above problems. First, the Landauer formula 
Q is a simple and clear description of the heat transport 
in purely ballistic regime at low temperatures. There 
are also a few other approaches @, HI with different ap- 
proximations such as the description based on density 
matrices. One of the features of the existing theories is 
that they work in either the ballistic regime or diffusive 
regime, but not both. It is rather difficult to have a 
complete theory that can encompass both regimes, apart 
from phenomenological treatments d, • 
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The nonequilibrium Green's function formalism de- 
scribed in this article is a serious attempt to be such 
a complete theory. This theory is certainly "first princi- 
ples" given the atomic potentials. Several similar formu- 
lations have already been done at the elastic level with- 
out nonlinear interactions [ll|, [IH, [TH, [l4j]. In Ref . (l5j . 
we presented a new formulation with the nonlinear in- 
teraction treated systematically. Mingo also gave a sim- 
ilar result [l6| for the nonlinear interactions. Our ap- 
proach is essentially a generalization of the nonequilib- 
rium Green's function method in electronic transport 
[l7l [Hj], with fermions replaced by bosons as basic en- 
tities. Although the techniques have been extensively 
described in the literature, most of them are centered 
around fermions. In this paper, we give a description 
of the method, emphasizing on both the theoretical for- 
mulation and computational implementation. We begin 
with an elementary introduction of the Green's functions, 
including the contour-ordered Green's functions. We 
then discuss equations of motion of the Green's functions, 
Feynman diagrammatic expansion, and Dyson equations. 
We discuss the heat current and derive a formula for effec- 
tive transmission when there are nonlinear interactions. 
After the introduction of the method, we present results 
of one-dimensional (ID) chains, benzene ring, and carbon 
nanotubes, and discuss some of the interesting features 
in such systems. 



II. NONEQUILIBRIUM GREEN'S FUNCTION 
METHOD 



A. Definition of Green's functions 

The nonequilibrium Green's function methods are dis- 
cussed in Refs. [HI EH, and equilibrium Green's functions 
are explained in many textbooks such as Refs. [13, [H|- 
For completeness, we introduce our notations and give 
a quick review of the definitions and properties of vari- 
ous Green's functions in this subsection. We define the 
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retarded Green's function as 

G r (t,i/) = -^e(t-t')([u(t),u(t') T })> (i) 

where the column vector u with component Uj is not 
the usual displacement operator from equilibrium but a 



renormalized one, i.e. 



such that the ki- 



netic energy is of the form (\/2)u T u. The superscript 
T stands for matrix transpose. The square brackets are 
the commutators. G r is a square matrix with elements 
G r jk (t,t>) = -(i/h)6(t - t'){[ Uj (t),u k (t')]). The phys- 
ical dimension of G r is time. By definition, G r (t,t') 
equals zero when t < t' . In equilibrium or nonequilib- 
rium steady states, the Green's function depends only 
on the difference in time, t — t' . The Fourier transform 
of G r (t - t') = G r (t, t') is defined as 



/+oo 
G r (t)e luJt dt 
-OO 



(2) 



Note that we use square brackets to delimit the argument 
for the Green's functions in frequency domain. The di- 
mension of the retarded Green's function in frequency 
domain is time squared. The inverse transform is given 
by 



G r (t) 



1 

2^ 



+ 00 



G r [uj]e- lut duj. 



(3) 



For notational simplicity, we'll set h = 1. We can always 
get the final required expressions by a simple dimension 
analysis. 

The rest of the definitions are the advanced Green's 
function 



G a (t,t') = i6(t'-t)([u(t)Mt'f]), 
the "greater than" Green's function 

G > {t 1 t') =-i{u(t)u{t') T ), 
the "less than" Green's function 

G < (t,t') = -i(u(t')u(t) T ) T , 
the time-ordered Green's function 

G'(M') =9(t-t')G > (t,t') + 9{t' -t)G<{t,t'), 
and anti-time-ordered Green's function 

G\t, t') = 9{t' - t)G> (t, t') + 6{t - t')G< (t, if). 



(4) 
(5) 
(6) 
(7) 
(8) 



The following linear relations hold both in frequency and 
time domains from the basic definitions: 



Gr s~ia f i > / 1 < 

G' + G* = G> +G<, 
G 1 - G l = G r + G a . 



(9) 
(10) 

(11) 



Also, the relations G r = G* - G< and G a — G < — G l 
are useful. Out of the six Green's functions, only three 
of them are linearly independent. However, in systems 
with time translational invariance, the functions G r and 
G a are Hermitian conjugate of one other: 



G>] = (G>])t. 



(12) 



So in general nonequilibrium steady-state situations, only 
two of them are independent. In this paper, we consider 
them to be G r and G < . There are other relations in the 
frequency domain as well: 

G < [w] t = -G<H (13) 
G r {-Lu] = G r [uj]\ (14) 
G<[-lu] = G > M T = -G < M*+G>] T -G r M*.(15) 

The last two equations show that we only need to com- 
pute the positive frequency part of the functions. 

In thermal equilibrium, there is an additional equation 
relating G r and G < : 



G<[uj] = f(cu)(G r [ij}-G a [ij]^, 



where 



i 



exp((3hui) — 1 



(16) 



(17) 



is the Bose-Einstein distribution function at temperature 
T = 1/(^8/3)- Equation P^|) is obtained by writing the 
Green's functions as a sum of energy eigenstates (known 
as the Lehmann representation). Thus in equilibrium, 
there is only one independent Green's function; we take 
it to be G r . 

The contour-ordered Green's function is a convenient 
book-keeping to treat the different Green's functions in 
a concise notation. We can consider a contour-ordered 
Green's function as a function G(r, r') with arguments r 
and t' defined on the complex plane. The contour runs 
from — oo slightly above the real axis to +oo and loops 
back from +oo slightly below the real axis to — oo. The 
contour-ordered Green's function can be mapped onto 
four different normal Green's functions by G aa (t,t') = 
lim £ ^ + G(t+iecr, t'+iecr'), where a — ±(1), and G ++ = 
G* is the time-ordered Green's function, G~~ = G f is 
the anti-time-ordered Green's function, G^ = G < , and 
G~ + = G>. 

In a noninteracting harmonic system with a Hamilto- 



Hn = —u T u H — u T Ku, 
2 2 



(18) 



the retarded Green's function in the frequency domain is 
given by 



G>] = ((io + itfl-K) 1 



(19) 



where / is an identity matrix, and rj — > + . Adding 
a small ry helps to choose correctly the inverse Fourier 
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transform integration path in the complex w-plane, so 
that the retarded Green's function has the required 
causal property, G r (t) = for t < 0. Equation (fP9|) can 
be derived by an equation of motion method. We note 
that the retarded Green's function is a symmetric matrix 
since K is symmetric, but this feature is not preserved 
with nonlinear interactions. 



B. Model and adiabatic switch-on 

Our system is a junction with a central region and 
two leads which serve as heat reservoirs. We treat the 
leads explicitly as quasi-one-dimensional periodic lat- 
tices. This setup is experimentally relevant and con- 
ceptually useful for computation. We consider non- 
conducting solid and treat only the vibrational degrees 
of freedom for heat transport. Let the (mass-normalized) 
displacement from some equilibrium position for the jth 
degree of freedom in the region a be uf; a — L,C,R, 
for the left, center, and right regions, respectively. The 
quantum Hamiltonian is given by 



H L +H c +H & +V+H m 



H=Y J H a + (u L ) T V LC u c + (u c fV CR u R + H n , 



(20) 



a=L,C,R 



where H a = h(u a ) J 



u a ) K a u a 1 u a is a column 



vector consisting of all the displacement variables in re- 
gion a, and u a is the corresponding conjugate momen- 
tum. K a is the spring constant matrix and V LC = 
(V CL ) T is the coupling matrix of the left lead to the cen- 
tral region; similarly for V CR . We note that the dynamic 
matrix of the full linear system is 



K 



K L yLC o 
yCL R C yCH 

V RC K R 



(21) 



The nonlinear part of the interaction will take the form 
H » = I E T ^ u ? u f u k + \ E T »kl «P- (22) 



ijk 



ijkl 



The quartic interaction is important in stablizing the sys- 
tem, as a purely cubic nonlinear interaction makes the 
energy unbounded from below. 

We need to answer an important question: what is the 
distribution of the system? The distribution enters the 
definition of the Green's functions as a density matrix 
p(0) for the average (•••). For an equilibrium problem, 
this is just the Boltzmann factor, but in a nonequilibrium 
situation, it is not known and must be computed in some 
way. The adiabatic switch-on gives us a clear concep- 
tual framework how this problem can be solved, at least 
formally. We imagine that at t = ™oo the system has 
three decoupled regions, each at separate temperatures, 
Tl, Tc, and Tr. The nonlinear interactions are turned 
off. The Green's functions g a are known and take the 
form of Eq. ([15]). The couplings V LC and V CR are then 




t = - 00 

Equilibrium at T t=T 

(=0 
Nonequilibrium 
steady state 

FIG. 1: A schematic to illustrate the two adiabatic switch- 
ons. 



turned on slowly, and a steady state of the linear system 
is established at some time T <C 0. The Green's functions 
of the linear nonequilibrium system will be denoted by 
Go- For this linear problem, the result does not depend 
on Tq- Finally, the nonlinear interaction H n is turned 
on, and at time t = 0, a nonequilibrium steady state is 
established, see Fig.Q]for illustration. The full nonlinear 
Green's functions will be denoted by G. 

The density matrices at time t = — oo, T, and t = are 
related in the following way in the interaction picture: 

p(T) = S (r,-oo)p(-oo)5 (-oo,r), (23) 

So (t,f) = Te- l ti V{t " )dt " ', (24) 
p(0) - S(0,T)p(T)S(T,0), (25) 

S(t,t') = Te -ip t , H ^") dt '\ (2 6 ) 
where T is the time-order operator. 



C. Equation of motion method 

In this subsection, we derive the equations of motion 
for Green's functions. The equations of motion can be 
used to perform systematic perturbation expansion for 
the Green's functions, or as a starting point for mean- 
field approximations. There are at least two ways to de- 
rive the equations of motion for nonequilibrium Green's 
functions. The first is to derive the equations of motion 
for the time-ordered Green's function and then gener- 
alize it to the contour-ordered version by evoking the 
structure isomorphism of the two sets of Green's func- 
tions [13]. Another possibility is to consider directly the 
contour-ordered Green's function [22]] . We'll take the lat- 
ter approach. 

We define a general n-point contour-ordered Green's 
function as 



G 



OL\ ,OL<2 . ■ ■ -,Ot r . 



(Vl,T2,- •• ,T n ) = 

-i(T T u^(n)u^(T 2 )...u^(r n )). (27) 



The index a = L,C,R labels the region, j labels the 
degrees of freedom in that region, and r is the contour 
variable. The function is symmetric with respect to si- 
multaneous permutations of the triplet (a, j, r). We view 
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the function as an analytic function (at least along the 
path) in variable t which varies on the Keldysh contour 
[23| . We also introduce a generalized 9(t,t') function 
which is defined to be 1 if r is later than r' along the 
contour and otherwise. Its derivative is a generalized 
5-function, 6(t — t') = a8 a ,a'5(t — t'). The 9 function can 
be used to represent the contour-order operator T T . For 
convenience, we allow for n — 0, which is just a constant 
—i. Our definition is slightly different from that of the 
usual quantum field theory where our prefactor is always 
—i (instead of the usual (— i)™/ 2 ). 

To start with, we consider one of the important cor- 
relation functions for heat transport, G c - L . The first 
derivative with respect to the second argument always 
leads to a direct differentiation inside the contour-order 
operator since the equal-time coordinates u commute, 

= -i(T T uf(T)^(r')). (28) 

The second derivative is similar in this case as u and u L 
also commute at equal time. Thus, substituting the equa- 
tion of motion for ul (having identical form in quantum 
and classical cases), we get 



d 2 Gf :i L (T,T<) 



-£ Gf^(r,r')K r L nl ^2 Gf£ (r,r')V r CL 



In matrix notation, it is 



(29) 



d 2 G c ' L {r, rQ 
dr' 2 



+ G c > l (t, t')K l = -G c > c {t, t')V cl . 



(30) 

To solve this differential equation on the contour, we de- 
fine the (contour-ordered) Green's function of the left 
lead to satisfy 



d 2 g L (T',T") 
dr' 2 



+ K l 9 l (t',t") = -I6(t'-t"). (31) 



We can think of — g L as the inverse of the operator 
d 2 /dr' 2 + K L . By multiplying 5 L (r', t") to Eq. (JSU]) from 
the right, G C ' L {t, t') to Eq. (JUl) from the left, and then 
subtracting the two equations and integrating over the 
variable r' along the contour, we get (after integration 
by part) 



G c > l (t,t") = / G C ' C (T,T / )V CL g L (T , ,T , ')dT' + 



dG C >\T,T') 

dr' 



9 l {t',t")+G c > l {t,t') 



,,dg L {r',T") 



- — OO — It 



- — oo-\-ie 



.(32) 



The reason that the terms in the second line are iden- 
tically zero involves some considerations. First, since 
the central part and left lead are decoupled at the time 
t = — oo + ie, we must have G C ' L (t, — oo + ie) = and 
-^ t G C ' L (t,~oo + ie) — 0. This takes care of the lower 
bound. Next, we require that the free left lead Green's 
function satisfies the boundary condition g L (r' , r") = 0, 



■gp-g L (r' ,t") = 0, at the upper limit r' — — oo — ie. 
This requirement is equivalent to <?*(— oo) = 0, and 
g > {— oo) = 0, or g r (— oo) = (as well as that their 
first derivatives are zero). The condition g r {~ oo) = is 
consistent with the causality requirement of the retarded 
Green's functions, but the conditions for g l or g > are in- 
sufficient, as the information about temperature of the 
system is unspecified. We fix the ambiguity by requiring 
that 



g>(t = 0) = -i(u L (0) U L (0) T ) 



(33) 



be the equilibrium value at inverse temperature (3l- In 
other words, g L should exactly be the contour-ordered 
Green's function of the free lead as defined in subsec. Hl Al 
Such a choice is consistent with the adiabatic switch-on. 
A full justification of the final result, 

G c > l (t,t")= f G C ' C (T,T')V CL g L (T',T")dT', (34) 



can also be given from a perturbation expansion, as has 
been done in Ref. [l7| . 

Equation (f3~4")) can be generalized to the case of G C,C,L 
with a similar result. We now consider the Green's func- 
tions of the form G ,C ' " ,C involving only the central 
part. The aim is to derive a recursion relation for the 
central part Green's functions by eliminating references 
to full Green's functions involving the leads. This type of 
Green's functions involves two new features, (1) we must 
take care of the derivatives of the 6 functions which pro- 
duce commutators [u c , ii c ), (2) we generate higher order 
Green's functions due to the nonlinear interactions. 

The contour order leads to nl terms, each of which is 
a permutation of the displacement variable u from the 
original order 1,2, ■■■,n, multiplied by a term of per- 
mutation of the arguments t from the standard order 
6{t\ — T2)9(t2 — T3) • • • 6{r n -i — t„). Differentiation of 
these 9 functions leads to a sum of contractions between 
the variable that we are differentiating and all the other 
variables. This reduces the order of the Green's function 
by two. We note that 



(35) 



For the moment, we'll consider only the cubic nonlin- 
earity. Four terms are produced when u c is substituted 
inside the order operator, due to linear coupling K , 
V CL , V CR , and nonlinear Tijk- The first three terms 
are similar in structure to that in Eq. (|29p . The cubic 
nonlinear interaction increases the value of n by 1. We 
note that the equation for u is evaluated at the same 
time t. This would be complicated if we have to keep 
track of which r's are equal in the Green's functions. We 
transfer this information to the coupling Ty^ to define 
the three-body interaction in contour variable form (for 
the force) 

jT l3k {Ty,T")^{T')u C k (T")dT'dT". (36) 

jk 
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FIG. 2: Recursive expansion rule for Green's functions. A 
vertex with n double lines denotes an n-point Green's func- 
tion. The single line denotes Go- A single-line three-terminal 
vertex is associated with T(t,t' ,t"). 

Since the contour integral is J dr — a f^°° dt, we 
must define 

Tijk(T,r' ,t") = T i jk(r'5 ( r t< r>6(t-t')a"5 ert(T '>6(t-t"). (37) 

This is because our (^-function has the standard property 
that J 5{t — r')dr = 1, but when written in the ordinary 
time t and index <r, we must have the extra a factor to 
get the required value of +1. 

The two terms involving the leads produce Green's 
functions of the form G ' ,L and G ' '"'' R , which can 
be replaced by G • using the result Eq. (f3"4"|) derived 
earlier. We can now solve for the center-only Green's 
function using the Green's function of the linear system: 

(£^+K c \ G (t,t') + J dr"{v CL g L (r, r")V LC + 

V CR g R (T,T")V RC )G (T",T>) = -I5(t-t<). (38) 

We then obtain the following simple recursive rules for 
Green's functions (as shown in Fig. [2]): 

(1) Replace leg 1 by inserting a nonlinear coupling 
T(t,t',t") such that the outer leg is immediately con- 
nected with Go while the other two terminals increase 
the order of the Green's function by 1. Quartic interac- 
tion is similar, but the process will increase the order by 
2. 

(2) Add imaginary unit i times a sum of the n — 1 
graphs formed by pairing each leg with leg 1 and con- 
necting with the propagator Go, multiplied by a (n — 2) 
order remaining Green's function. 

(3) Symmetrize the graphs, if desired, i.e., do steps (1) 
and (2) for every leg, 1, 2, • ■ •, n, add them up, and then 
divide by n. 

The above rules may be conveniently implemented in 
a symbolic computer language, such as Mathematica. 

D. Feynman diagrams and Dyson equations 

The contour-ordered Green's function can also be ob- 
tained by a perturbation expansion of the interaction pic- 



ture evolution operator (scattering matrix operator) and 
is expressed as: 

G jk (r,r') = -i{T T uf{T)u%{r')) 

= - l {T T u]{T)u I k {r')e~^<^'')^ (39) 

where the displacements refer to the central region; we 
have dropped the superscript G for brevity. The opera- 
tors in the top line are in the Heisenberg picture, while 
that in the second line are in the interaction picture. The 
Green's function Go of the linear system (when H n = 0) 
can be computed from that of the free subsystems (an 
integral equation form of Eq. (j3"8|) ): 

G {t,t')=9 C {t,t') + JdT 1 dT 2 g C (T,T 1 )J:(T 1 ,T 2 )G (T 2 ,T'), 

where E = E^ + E#, 

E L = V CL g L V LC , (41) 

similarly for E#. This Dyson equation can be derived by 
considering the first step of the adiabatic switching-on 
process. Since it is a linear system, the self-energy E is 
known exactly. 

By expanding the exponential in Eq. (139|) . a series in 
the nonlinear interaction strength is obtained. The re- 
duction in terms of the unperturbed Green's function Go 
relies on the fact that Wick's theorem [24[ is applicable 
here. Although a formal proof of Wick's theorem is dif- 
ficult, we can see that the density matrix p(0) must be 
quadratic in exponential. This is because the p{— oo) is 
quadratic and the system is linear. The fact that the 
equation of motion method and the perturbation expan- 
sion give identical results is a confirmation of the validity 
of Wick's theorem. 

The expansion contains connected as well as discon- 
nected Feynman diagrams. The disconnected diagrams 
are constant in time, and give rise to a thermal expan- 
sion effect. We can show that these diagrams do not 
contribute to the thermal transport, as they are propor- 
tional to S(oj) in the frequency domain. The thermal 
current formula has a factor of uj which makes it zero. 

Finally, the connected part of the Green's function sat- 
isfies a similar contour-ordered Dyson equation relating 
G c to Go through a nonlinear self-energy E n : 

G c (t,t') = G (t,t')+[> in dr 2 Go (r,n ) E„ ( n ,r 2 ) G c (r 2 ,r' ) . 

(42) 

In ordinary Green's functions and in frequency domain 
(u> argument suppressed), the Dyson equations have so- 
lutions 

Gl = Gg f = ((cu + ttfl-KC-^y 1 > ( 43 ) 
G< - GSE<Gg, (44) 

G r c = (gS^-E;)" 1 , (45) 
G< = G£(E< + E<)G«. (46) 



6 



O 



+ 2 i 



9 



+ (-! 




we have 



3G 



33 ,34 ,35 -J6-XT4 



+ (-4) V' + (4) V + C-2)^V^ + 3 



(-6) 



OO + (12) "0- + ( ~ 12) .rv <- 6 i(^i + (- 6 > 



+ (-12)- 




- (-6) , 



+ (-6) ( V (-6) 



6 *»qp 



<«JS.«'-«)ciS3,«»»4** ' 49 ' 



The result can now be transformed to frequency domain 
if desired. 



E. Thermal current and conductance 

We define the current flow from the left lead to the 
central region as 



FIG. 3: The Feynman diagrams for nonlinear self-energy E„ 
with the prefactors for graphs of cubic and quartic interac- 
tions to second order in Ti. 



Il = -(H L (t)). 



(50) 



The connected part of the Green's function G c will be 
used for G in Eq. (|52[) below, and the extra subscript c 
will be dropped from this point. 

In Fig. [3] we give the Feynman diagrams to "second 
order" of the nonlinear interactions for cubic and quartic 
terms, where the graphs have explicitly H and H 2 in the 
coefficients; the neglected graphs have higher powers of %. 
These diagrams are the same as those in Ref . [25| , but the 
interpretation is different. It is for the contour-ordered 
Green's functions here and for the retarded Green's func- 
tions of the creation/annihilation operators in Ref. [2S| . 
To find the actual expressions for each of the diagrams, 
we have the following simple rules: (1) label each interac- 
tion vetex line by the r variables and site indices, with the ■ n 
interaction function Ty fe (rj,Tj-,T fc ) (or Tijkl(Ti,Tj,T k ,Ti) G bcM V 9l M 



The interpretation of the current is somewhat subtle as 
the leads are semi- infinite in extent. The definition is 
meaningful. By the Heisenberg equation of motion, we 
obtain, at t = 0, I L = ((u L ) T V LC u c ). The expecta- 
tion value can be expressed in terms of a Green's func- 
tion G^ L (t,t') = -i(u L (t')u c (t) T ) T . Using the fact 
that operators u and u are related in Fourier space as 
u[lo] = {—iuj)u[uj\, we can eliminate the derivative and 
get, 



LC 



;]) u) duj. 



(51) 



Using the result derived earlier relating the mixed lead- 
center Green's function to the center-only Green's func- 
tion and applying the Langreth theorem [HI, [26[ to write 
contour-ordered Green's function in ordinary Green's 
functions in frequency domain, we have Gp L [w] = 

]V CL g a L [u\. The final expres- 



G'cc 



to 



for the quartic interaction). (2) Each line connecting the 
labels is associated with a propagator Go- (3) Integrate 
over all the internal variables along the contour; sum 
over the internal site indices. Most of the integrations 
are easy, as Tijk(- ■ •) contains two ^-functions in r. For 
example, the second graph is 



sion for the energy current is 



i r + °° / 

I L = / duuiTrl G r [oj]T,f [oj] + G < [uj] E£ [u>] 

27T/-00 V 



(52) 



E 



dT 3 dT 4 dT 5 d,T 6 T juj2tj3 (n, r 2 , r 3 ) 



33 ,34 ,35 ,3e 

G o ,3334, ( T 3 , Ti)T u j 5 j e (r 4 , t 5 , t 6 ) Go j 5je (t 5 , t 6 ) , (47) 

where t\ and t 2 are external variables and the rest are 
internal. After integrations, we get 

E / dr i T 3i ,32 ,j3<5( T 1 _r 2 )Go ,33 34.i T l , T i)Tj 4 ,j 5 Jfpojsje ( T 4 , T 4 ) 

,34 ,35 ,36 

(48) 

Now we can rewrite in real time t and the branch in- 
dex a using the mapping G(t, t') — > G CT<T (t,t'), J dr — > 
E ff =±i X^o °" d *> and S ( T - T ') —> vS a ,cr'S(t - t'). Noting 
that Gq(t4,T4) = Gq(0) by time translational invariance, 



where the self-energy due to the interaction with the lead 
is Yil = V CL giy LC . For simplicity, we have dropped the 
subscript C on the Green's functions denoting the central 
region. This result is the phonon analog of the electronic 
current formula [Tt], HH • 

Since the energy current must be conserved, we have 
Il = —Ir. In addition, Il is real. We can obtain a 
symmetrized expression for the current 



1 

~ 4 
lo Tr 



{(G 



II- Ir- Ir) = 



1 



4tt ,y 

G a )(£<-£<)+zG<(r fl -r L 



did 



where r a = i(E£ - E£). 

We define the thermal conductance as 

- r 1 61 



)}, (53) 



(54) 
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FIG. 4: Mean-field approximation to the self-energy. The 
double line denotes the full Green's function G; the single 
line is for Go- 



where AT is the difference of the temperatures in leads, 
such that T L = T + AT/2 and T R = T- AT/2. Using 
Eq. ([55)1 . we can actually take the limit AT — > by in- 
troducing variational derivatives of the Green's functions, 
self-energies, or current as, e.g., 



SG 

— — = hm 

5T AT^O 



G(T L ,T fi ) - G(T,T) 

Tr - Tr 



(55) 



The difference in the numerator is the value of a func- 
tion when the leads are at two different temperatures Tl 
and Tr, respectively, minus the value when the system 
is in equilibrium at Tl = Tr = T. With this notation, 
we can express the conductance in a form similar to the 
Landauer formula [27] for the ballistic transport as 



1 

27 



(1ujujT[uj\ dT 



(56) 



with an effective transmission coefficient, 



f[u 



l -Tr{G r {T L + l -T r , 



S)G a T R ) 



l -Ti{G a T L G r {T R + l -T n + S)}, 



(57) 



where the nonlinear effect is reflected in the extra terms, 
r„ =i(E;-E»), and 



s = % 



/( 



6T 



ST ■ 



ST 



dT 



(58) 



Equation (|57|) is a generalization of the Caroli formula 
[281 ] for ballistic transport. We'll discuss in Sec. IIII Bl 
how the variational derivatives of the self-energy can be 
computed. 



F. Mean-field approximation 

We have seen in Ref. [151 ] that high-order Feynman di- 
agrams are important at high temperatures. However, it 
is computationally very expensive to include these high- 
order diagrams, although we can do calculations on very 
small systems such as a three-site model. The standard 
technique is to use a mean-field approximation which in 
some sense partially takes into account these high-order 
graphs. We propose an approximation to the self-energy 



shown in Fig. [4] where the double line represents the non- 
linear full Green's function G, while a single line is the 
usual bare Green's function Go- With such a choice, we 
find that to next order, all the cubic interaction graphs 
are reproduced exactly (including the prefactor), except 
that two of the diagrams are missing. 

With the self-consistent mean-field approximation, the 
solution has to be found iteratively. The convergence of 
the iterative process becomes an issue. With only cubic 
interactions, the system of equations will not converge at 
high temperatures. We understand that this is simply 
a consequence of the unstable nature when only cubic 
terms are used. We thus add the fourth order potential 
which is essential for stability. However, the convergence 
is still difficult at high temperatures, mainly due to the 
singular nature of the Bose distribution. 



G. Classical limit, molecular dynamics 

At high temperatures, the classical molecular dynam- 
ics (MD) is a correct and numerically exact method for 
the heat conduction problem. However, the standard 
heat baths used in MD such as Nose-Hoover do not cor- 
respond to exactly the action of the leads in our models. 
A correct heat bath should follow a Langevin dynamics 
with colored noises. Dhar [291 ] has given a derivation for 
the junction problem. The idea is to solve the lead de- 
grees of freedom u L and u R formally first, which is easy 
since it is a linear system. The result is then substituted 
into the equation for the center region. This gives the 
following equation: 

u c = -K c u c + F n (u°)-I ' ^ r {t,t')u c {t')dt' + Ct + 

Jt 

(59) 

where F n is the nonlinear force, E r is the same retarded 
self-energy of the lead, XT = TT L + Y7 R , used in the 
nonequilibrium Green's function calculation, but in the 
time domain. By integration by part, we could also cast 
Eq. (|59p in a standard generalized Langevin form (30| 
involving a damping force, but K c has to be shifted to 
K +T, r [uj = 0]. An extra contribution from the left lead 
is 



tL(t) = V CL (g L (t,T)u L (T) 



(60) 

where g L is the retarded Green's function of the left lead 
in the time domain. The expression for the right lead 
£_r is similar. A time T in the lower limit in the inte- 
gral and Eq. (|60p is some time in the remote past where 
the system is decoupled and the leads are in respective 
thermal equilibrium. We'll take the limit T — > — oo. We 
turn Eq. (|59|) into a stochastic differential equation by 
requiring that u L (T) and ii L (T) are random variables 
satisfying the Boltzmann-Gibbs distribution 



p(u L ,u L ) oc e-MK^r^+M^r*^). 



(61) 
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Thus, £l is colored noise with the properties 

(&(*)> = 0, (62) 
/?L(a(t)a(t') T > = V CL (gl<t,T)jLgl(?,T) T 

+ g r L (t,T)g r L (t\Tf)v LC (63) 

O 

X r L (t",t')dt", (t>t'). (64) 



The last line of the above equation is a kind of 
"fluctuation-dissipation" theorem. We note that the dot 
in g denotes the derivative with respect to T. For a sen- 
sible heat bath, the noise should not depend on T, but 
only on the difference t — t' . Indeed, this can be done if 
we rewrite the result in vibrational cigen-modes. We can 
write the Fourier transform of the noise correlation as 



Im2V 



CL 



(K L ) 



-1/2. 



f 



irj)I - K L 



V 



LC 



(65) 



where the square root of a matrix is defined by its eigen- 
value decomposition taking only the positive eigen fre- 
quencies. This is well defined since K L is positive defi- 
nite. The properties of the right lead noise £r are anal- 
ogous. Such a set of stochastic differential equations can 
be simulated on a computer. 



III. IMPLEMENTATION DETAILS 
A. Surface Green's functions 

The Green's functions of the free leads g r are required 
in calculation for the lead self energies S. Although the 
lead Green's function is a semi-infinite matrix, due to a 
localized coupling V or V CR , only a corner set of ele- 
ments is required. These submatrix elements are termed 
surface Green's functions. The surface Green's function 
can be computed rather quickly by recursive iterations 
[3~fl ] or decimation. For completeness, we give an algo- 
rithm in pseudo-code form. The reader is referred to the 
literature for the derivation of such algorithm. 

The program takes three square matrices of equal size 
for fcoo, fcn j an d fcoi- The matrix fcoo is the block imme- 
diately adjacent to the center, while fen, and fc 01 = kf a 
are repeated for the semi-infinite chain of the lead. They 
form the whole dynamic matrix of the lead in the follow- 
ing way: 



K 1 



The left arrow 
error tolerance. 



/ fcoo fcoi • • • 
fcio fcn fcoi 
fcio fcn fc i 



(66) 



/ 



V ■ • ■ fcio 
below denotes assignment; e is an 



S <r- 


- fcoo 


e <- 


■ fcn 


a <- 


- fcoi 




-a T 


g «- 


-[(" 


do 





■vifl-e]- x 
agP 



s *— s - 

if \s' — s\ < e exit 
e' <— e + ag/3 + (3ga 
a' <— aga 
(3' (3gf3 
s «- s' 
e^e' 
a <— a' 
(3^0' 

g 4- [(u + ir)) 2 I - e] 
end do 

g r <- [{u + irffl- s]- 1 . 



B. Fast Fourier transform and variational 
derivatives 

The diagrams in the mean-field approximation can be 
computed in the time domain easily and efficiently; the 
convolutions in frequency domain become simple mul- 
tiplications in time domain. Thus we have adapted a 
fast Fourier transform method for the calculation. This 
reduces the computational complexity, particularly for 
the quartic interaction terms. Typically 10 3 Fourier grid 
points are needed for good convergence. Using the Mes- 
sage Passing Interface (MPI), a parallel implementation 
with data partition is made. Good linear speed-up is 
obtained upto 64 processes. 

In the mean-field treatment of the heat transport, there 
is a need to find the variational derivatives of the Green's 
functions and self-energy for the conductance calculation. 
This is done by solving iteratively the set of equations 



5G r = G r 6E r n G r , 
5G< = G r 5S:G < - 



G r (<5£ < +<5£<)G a - 



(67) 

-G<c5££G a (68) 



The expressions for iJE^ and 5S< can be obtained by 
differentiating the self-energy results directly. For con- 
vergence control, linear mixing of the new and old values 
is used. 



IV. APPLICATIONS 

A. One-dimensional cubic onsite model 

One important question to ask about the nonequilib- 
rium formulation is that whether the theory gives bal- 
listic and diffusive heat transport in a single framework. 
Clearly, if the nonlinear self-energy £„ can be computed 
accurately, there is no doubt that the answer is yes. How- 
ever, this is not clear if only the leading order perturba- 
tion theory is used. In Ref. [32j , Spohn gives a derivation 
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of a Boltzmann equation for the phonons in a simple ID 
chain with onsite cubic nonlinearity. It is clear that in 
the Boltzmann equation framework, we obtain Fourier's 
law and diffusive heat transport. We demonstrate in this 
subsection that diffusive heat transport is not observed 
if only perturbative result is used. This subsection also 
serves as an example of illustration of the concepts dis- 
cussed earlier by a simple and concrete model. 

We consider a ID chain with inter-particle spring con- 
stant K and onsite spring constant K . We take a sim- 
ple cubic onsite potential for the interaction, gi^Zjti^, 
for j belonging to the central region. This means that 
Tijk = tdijSik- Apart from the nonlinear interaction, the 
left lead, central region, and the right lead are identical. 
The classical equation of motion is 



Kuj-! + (-2K - K )uj + Ku 



3+1 



t in 



(69) 



where tj = t for 1 < j < N and tj = for j < 1 or 
j>N. 

The retarded Green's function Gg can be obtained by 
solving the linear equation + irj) 2 — K]Gq = I, where 
(infinite in both directions) matrix K is 2K + Kq on the 
diagonal and —K on the first off-diagonals. The explicit 
solution is 



Gr 
Ojl 



\3-l\ 



(70) 



(Ai-A 2 )A-' 

where Ai and A2 are the roots of the equation 

KX- 1 + fl + KX = 0, fl = {u; + irf) 2 -2K- K . (71) 

Note that A1A2 = 1. We define the two roots such that 
|Ai| < 1 and |A 2 | > 1. 

The surface Green's function g r satisfies a similar equa- 
tion as Gg except that it is semi-infinite in extent. We 
consider the left lead (j < 0). The result for the right 
lead is identical. Since the matrix V LC is nonzero only 
for one corner element, we need the go = goo component 
of the matrix g r . Consider only the j — column, the 
equations for g r in component form are 



figo + Kg-i = 1, 



(72) 



Kg j - 1 +ilg j +Kg j+1 = 0, J = -1,-2,- • • (73) 
Substituting the trial solution gj = cX^ , we find that 



Ai 

90 = -K- 



(74) 



The "less than" Green's function for the left lead is then 
(c.f. Eq. ((Til)) 



9l = 2i/ilmg = -2if L 



Im Ai 

K ' 



(75) 



where fi is the Bose-Einstein distribution function at 
the temperature of the left lead. The Gq Green's func- 
tion can be obtained with the equation Gq = GqS < Gq, 



E = S L + S L = V gjJV , similarly for E fl . 
The matrix elements of E < are all zero except that 
En = K 2 g^ and Ejvat = K 2 g^. Using these results, 
after some calculation, we can write 



GofM 



hxt 1 + f R X[~ j 
(Ai - X 2 )K 



(76) 



The Q(u>) function is defined to be 1 inside the phonon 
band and zero outside the band; more precisely, Q(u>) = 1 
if Kq < <jJ 2 < 4K + Ko, and is otherwise. 

Due to the onsite cubic nonlinearity, the leading or- 
der nonlinear self-energy calculation is simple. The first 
graph contribution is 



(77) 



The contribution from the second graph is a constant in 
2 £„?f'M = 2it 2 a5 a ^8 jl ct"G ^"[0]G :; ct "(0). 

a" .s 

(78) 

Both the summation over a" and over index s can be 
performed analytically. We find 



2 E„^ [u] = 2tt 2 <jS a ^S jl G r J [0]G t (0), (79) 



where 



G*M = 



G*(0) 



1 + Ai - X{ Af 
(l-Ai)(Ai-A 2 )A- ' 

l + (/L+/fl)6H 

{X l ~X 2 )K ' 

2^ 



(80) 
(81) 
(82) 



For numerical accuracy, we compute the conductance 
directly instead of using finite differences. Let the tem- 
perature on the left lead be T L = T + AT/ 2 and right 
lead be T R = T - AT/2. We define the thermal conduc- 
tance as that in Eq. (|54|) and the effective transmission 
by Eq. f|56[) . For the expression for effective transmis- 
sion T[u>], we also need to take the limit AT — > 0. The 
advantage of doing this is that the leading terms corre- 
sponding to equilibrium Green's functions cancel exactly, 
since there should be no heat current when Tl = Tr. 
We can also equivalently define the temperature varia- 
tion through the expansion: 

0)G r 

G r = G r cq + — AT + 0(AT 2 ), (83) 

and similarly for G < . We can simplify the effective trans- 
mission as 



f[u] = iTr{(G'--G )(S< 



iG < (T 1 



)}/(/! 



f R ) 



iK{Im\ 1 ){FL 1 -F* N ), 



(84) 
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FIG. 5: Thermal conductance of the cubic onsite model as a FIG. 6: Thermal conductance as a function of temperature 

function of temperature for several values of the cubic cou- for several system lengths N with t = 0.5 eV /(A 3 u 3//2 ). 
pling strength t = 0, 0.2, 0.5, 0.7, 1.0, and 2.0 eV/(A 3 u 3/2 ). 
The length of center chain is TV = 5. 



where the F function is defined as 



F = -(G r - G a ) 



S(G r - G a ) SG* 



ST 



5T 



2i 

dT 



F R is the same except that the first term is negative. In 
evaluating the F functions, we should set Tl = Tr = 
T as we have already taken the AT — * limit. The 
variations required in Eq. (|85[) can be computed by the 
variation in nonlinear self-energy, which in turn can be 
computed by the variation in the linear Green's function 
Go. Due to the Dyson equation, the variation of the 
retarded Green's function has a simple result: 



The variation in G < is rather complicated: 
SG< 



(86) 



SG < 



ST 



-AT = G r {E;<5G< 



6Z r n [G< + G r (X<G a + E;G<)] } - h.c. 
+ SG< + G r (<SE< + E^GtfE-)G°, (87) 

where Gq = G < (/-|-S°G a ), h.c. stands for the Hermitian 
conjugate of the preceding term, and SGq = 0, 

ST 2{\ X -\ 2 )K K 'dT' 



5E n ^ < [w] 



4i t 



2 r+ao 



du , G ^ < [u-u/]6G f l [u'].{89) 



of the thermal conductance R of a short nonlinear chain 
of length N — 5. The t — curve is the ballistic Lan- 
dauer formula result, the rest of the curves show the effect 
of nonlinearity. The prominent feature here is that the 
conductance decreases quickly as the nonlinear coupling 
t increases from to 2. In particular, the system behaves 
like a thermal insulator if the on-site nonlinearity is large 
enough. The thermal conductance is small at low tem- 
peratures and at high temperatures, as was found exper- 
imentally. At the extremely high temperatures, the con- 
ductance decreases with temperature according to 1/T 
(data not shown). This agrees with several other theo- 
ries of thermal conductance at high temperatures. 

Can we see diffusive transport in this model (under 
the perturbative approximation)? The answer is no. In 
Fig. [HI we show the size dependence for a fixed t. At tem- 
peratures below 200 K, the results are essentially inde- 
pendent of the lengths, showing ballistic heat transport. 
At around 2000 K, from sizes 2 to 10, we see decrease of 
the thermal conductance roughly as 1/N, but this trend 
does not continue, and the conductance saturates to some 
value for systems larger than 20. This feature is generic; 
the FPU model with first-order perturbation treatment 
also has a similar behavior. The curves also display some 
oscillatory structure. The number of oscillations appears 
to be roughly equal to half of the chain length N. The 
feature is damped out as the chain becomes longer. It 
would be nice if a molecular dynamics simulation could 
confirm these features. 



We have made analytic calculations as far as pos- 
sible. The rest of integrations and matrix inversions 
are calculated numerically. We use physical units with 
K = 0.625 eV/(A 2 u), K = Q.IK, (u for unified atomic 
mass units). Figure [5] shows the temperature dependence 



The failure to see diffusive behavior may be under- 
standable as our result takes into account only the first- 
order perturbation in the nonlinear interaction. Other 
features seem qualitatively correct in comparison with 
real systems; this is very encouraging. 
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FIG. 7: Thermal conductance of a Fermi-Pasta-Ulam chain 
with chain lengths 4, 8, 16, 32, 64, 128, 256, from top to 
bottom. Semi-infinite one-dimensional harmonic chains serve 
as heat baths with spring constants fc_L = 2.31 and ku = 2.25 
eV/(uA 2 ) for the left and right leads, respectively. Small 
onsite quadratic potentials are applied with fc£ nsltc = 0.006 
and fc?T ito = 0.01 eV/(uA 2 ). 



B. Fermi-Pasta-Ulam model 

The Fermi-Pasta-Ulam (FPU) model has been studied 
as a prototype model in nonlinear dynamics [331 ]. In re- 
cent years, a lot of works have been done on the thermal 
transport in the FPU model [H, (3|| . By molecular dy- 
namics and mode-coupling theory, Lepri et al. (3ol [37j 
found that the thermal transport in such a system does 
not have an ordinary diffusive behavior in the sense that 
the thermal conductivity diverges with the system size 
as some power, L a . The actual value of the exponent a 
is somewhat controversial. The earlier studies gave 2/5, 
but the latest analysis suggests 1 /3 [H, . 

Quantum thermal transport in FPU model has seldom 
been studied. In particular, the change from ballistic to 
diffusive transport has not been addressed. This ques- 
tion can not be correctly answered within molecular dy- 
namics approach as the ballistic transport is quantum- 
mechanical in nature. 

We use realistic model parameters derived from ex- 
panding the Morse potential, V{x) = D(e~ ax — l) with 
D = 3.8 eV, a — 1.88 A -1 , up to the fourth order in 
x = (ui — Ui + \) / ^Jm, where Ui is mass-renormalized dis- 
placement for the «th particle. It is a general FPU model 
with both cubic and quartic interactions. The leads do 
not have nonlinear interactions. We use the mass of car- 
bon, m = 12 u, in the calculation. Thus we can consider 
the results correspond to that of a carbon chain. 

In Fig. [7J we present the mean-field results of con- 
ductance function of temperature for different 
lengths of the chain. As we can see, below about 300 K, 
the conductance is independent of chain lengths; this is 
the ballistic region of heat conduction. At high tem- 




T (K) T (K) 

FIG. 8: Molecular dynamics results of the thermal conduc- 
tance for the model as in Fig. [7J (a) "Generalized Langevin" 
thermal heat baths from sizes 4, 8, 16, ... , 2048, from top to 
bottom, (b) Nose-Hoover thermostat with size 8, 16, 32, . . . , 
1024. 



peratures, the conductance decreases with the system 
size. For a completely diffusive heat transport, we should 
have R oc 1/L. We are far from this. Since in FPU 
model real diffusive behavior is impossible, the reduc- 
tion in conductance in the classical regime should be 
L Q_1 L~ 6 . The observed exponent for sizes up to 
256 is much smaller. In fact, this result is expected be- 
cause the mean free path in a typical carbon systems 
(such as carbon nanotubes [4(| or diamond) is of order 
fim. Since the conductivity decreases with temperature 
typically as 1/T, even at few thousand kelvin, the mean 
free path is expected to be 100 to 1000 lattice spacings. 
Thus we see only tendency towards diffusive behavior. It 
is also clear that ballistic transport and diffusive trans- 
port are at quite different time and length scales. The 
present calculation is already very large on an IBM 500- 
processor parallel computer. It is out of the possibility 
to see diffusive behavior for large sizes within the present 
approach. But this example does give us an indication 
that the mean-field theory could give diffusive behavior 
for appropriate systems. 

It is instructive to compare the mean-field results with 
molecular dynamics results. In Fig. [5J we present the 
MD results with the generalized Langevin-like heat-bath 
and the Nose-Hoover heat bath. The colored noise is 
implemented using a standard spectrum method with 
fast Fourier transform [41]. It turns out that the re- 
sults for the thermal conductance is very sensitive to 
the bath used and detail couplings between the bath 
and system. It is clear that quantum nature shows up 
already as high as 1000 K for this model chain. This 
can be understood by a very high Debye temperature 
T D = huj mllx /kB ~ 2000 K, where w max = \J^K C . The 
Landauer formula, Eq. (|56|) . provides a rigorous upper 
bound for the thermal conductance (when T[uj] = 1, 
T — > oo), K max = ^ B w max « 0.65 nW/K. Interestingly, 
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FIG. 9: Effective transmission of the benzene ring junction at 
a temperature of 300 K. The solid line is for a full model with 
force constants calculated from Gaussian, while the dotted 
line is from a simplified 2D model. The vertical bars indicate 
the vibrational frequencies of isolated benzene rings of the full 
model (topmost) and simplified model (dashed), respectively. 

this is satisfied at low temperatures by the generalized 
Langevin data. 

To what extent the mean-field theory is "correct" ? We 
can see that it is qualitatively correct for the accessi- 
ble temperature and size ranges. However, if we trust 
the Langevin heat-bath result, then, it appears that the 
mean-field approximation somehow over estimates the ef- 
fect of nonlinearity. 



C. Benzene ring 

Molecular devices have attracted a lot of attentions in 
recent years, due to their potential for future generation 
electronics. Experiments have been conducted on such 
systems; theoretical analyses of their electronic and ther- 
mal transport properties have been performed 0, |42j ■ 

In this subsection, we examine the thermal trans- 
port properties of a benzene molecule connected to two 
polyethylene leads. Two opposite sites of a benzene 
ring, which are normally the hydrogen atoms, are now 
bound covalently to carbon atoms, forming a 1,4-dibutyl- 
benzene molecule representing the whole system. The 
system has 14 carbon atoms and 22 hydrogen atoms. The 
interactions are calculated using Gaussian 03 [43j]. We 
use the Hartree-Fock method with the 6-31G basis set. 
We have also used a convenient new feature of Gaussian, 
which can calculate anharmonic interactions up to the 
4th order. 

Fig- HI (solid line) shows the effective transmission func- 
tion T[ui] of the benzene junction model at a tempera- 
ture of 300 K calculated using the mean-field theory. The 
peak features are complicated, resulting from the vibra- 
tional modes of the molecule, the interaction with the 



leads, and the vibrational spectra of the leads, as well 
as the temperature effect. To understand the transmis- 
sion function better, we correlate the peak features with 
the vibrational spectrum of an isolated benzene molecule 
calculated from Gaussian, shown in Fig. [9] as solid bars. 
The transmission peaks are roughly located at where the 
vibrational frequencies appear, but somewhat shifted to 
low frequencies. It is difficult to identify the peaks with 
eigenmodes unambiguously. 

We notice that the hydrogen atoms do not play much 
role in the thermal transport. In particular, there are 
frequencies around 3000 cm -1 associated with the vibra- 
tions of H atoms, but there is no transmission in these 
high frequencies. Vibrations of H atoms have only local 
effects. We confirm that the H atoms are less important 
by a much simplified model retaining only the degrees 
of freedom of the carbon atoms in a 2D model. The 
leads are pseudo polyethylene chains; we use an united 
atom approximation, so the CH2 group is treated as a 
single atom. The force constant between each atom is 
obtained from the DREIDING model 0. There are 
only bond interactions in the leads, and we have used the 
Morse function to calculate the bond potential, which has 
the form £>( e -"(«~«o) _ 1 )2 ) where D = 70 (kcal/mol), 

a = y/5 A , and R — Rq is the displacement (with the 
equilibrium bond length Ro = 1.4 A). For the benzene 
ring, we also include an angular potential, which has the 
form !(cos(0y fe ) + \)E, where E = 100 (kcal/mol). We 
make a Taylor expansion of the potential and calculate 
the force constants up to the 4th order. 

The transmission for the simplified model (dotted line 
in Fig. [9|) is smoother due to much reduced degrees of 
freedom. It is qualitatively in agreement with the full 
model. In this case, the peaks appear to be in good agree- 
ment with the vibrational modes of a hexagonal ring. 



D. Carbon nanotube junction 

Thermal transport of carbon nanotubes has been stud- 
ied by molecular dynamics lip . I4H |47| . |48| , Landauer 
formula ballistic transport }l3j |. and a phenomenologi- 
cal theory taking into account the nonlinear effect [49(. 
There are also experimental works on carbon nanotubes 
[50l ]. We have also studied the heat transport proper- 
ties of the (5, 0) carbon nanotube using a tight-binding 
(TB) method [51(, of which the TB parameters for the 
carbon atom have been checked recently [52^ . The left 
lead, the central junction, and the right lead consist of 
forty, twenty, and forty atoms, respectively. A unit cell 
of the (5, 0) nanotube has twenty atoms. The tube's 
axial dimension is long enough (about 20 A) to allow 
a single T point to be used for the /c-point sampling. 
We have written a force subroutine in order to perform 
atomic relaxation of the nanotube as well as to calcu- 
late force constants by numerical finite differences. The 
spring constants K a , the nonlinear force constants T^-fe 
and Tijki are systematically obtained by displacing one, 



13 



3 



5 2 
c 







1 

- 


1 1 1 




ballistic 


- 


nonlinear 




1,1, 



1^ ' 1 ' 1 ' ' — 

200 400 600 



T(K) 

FIG. 10: Thermal conductance function of temperature 

T for a twenty-atom (5,0) carbon nanotube junction. 

two, and three atoms, respectively. For each displace- 
ment, the forces acting on the atoms are used in the 
central-difference scheme to calculate the force constants. 
Since the TB method we have adopted is free from nu- 
merical noises, we have used a reasonably small maxi- 
mum displacement of /i ma x of 3 x 10 -4 A (i.e., the atoms 
are displaced not more than /i max from their equilibrium 
positions). We note that this value of h max is two or- 
ders of magnitude smaller than that commonly used in a 
supercell force-constant method [EH, (54[ ■ 

In Fig. [10] we present the results of thermal conduc- 
tance calculations. We find that for such a complicated 
system, the mean-field iterations are hardly converging, 
thus we only give the ballistic (where the nonlinear inter- 
actions are turned off), and a first-order perturbative re- 
sult using the (single-line) diagrams in Fig. [4j It is clear 
that the nonlinear interactions have a large effect even 
at moderately high temperatures. The thermal conduc- 
tance is slightly peaked around 600 K for the nonlinear 
result. We expect that the result should be reliable at 
least up to room temperatures. 

V. CONCLUSION 

We presented a nonequilibrium Green's function 
method for thermal transport in nano-junctions. We 
gave a unique derivation of the equations of motion of 
Green's functions and rules for generating and comput- 
ing Feynman diagrams in a contour-ordered Green's func- 
tion setting. The nonlinear interaction was treated with 



a systematic perturbative expansion for self-energies. For 
practical calculations, the mean-field approximation has 
to be adopted. We are still left with a technical dif- 
ficulty that the mean-field equations sometimes do not 
converge. Several versions of computer programs have 
been written for consistency checks (a Mathematica note- 
book for an initial fast prototyping, a Fortran program 
with MPI for very general models, and a specialized pro- 
gram for the ID cubic onsite model). We have applied 
the programs to a number of model systems. Several 
sophisticated 'first-principles' methods were used in the 
force-constant calculations for realistic models. For the 
FPU model, we can accurately predict the ballistic ther- 
mal transport and crossover to diffusive behavior at suf- 
ficiently high temperatures. We proposed to use gener- 
alized Langevin heat baths, which are more appropriate 
for nano junctions. For the benzene ring model, we ob- 
served interesting behaviors in the transmission function 
related to the atomic vibrational modes of the molecule. 
For carbon nanotubes, it is clear that nonlinear effect is 
quite strong even at low temperatures. 

The Green's function discussed here is more funda- 
mental than the phonon distribution function in the 
Boltzmann equation. The connection of the present for- 
mulation with the quantum Boltzmann equation or the 
phonon Boltzmann equation is something worth pursu- 
ing. Formally, they are related by the Wigner transfor- 
mation. If we can reduce the matrix of the Green's func- 
tions to vector of distribution functions (in momentum 
space) , potentially much larger systems can be treated if 
we are willing to forsake some atomic details. Clearly, the 
nonequilibrium Green's function method is at its best for 
small nanoscale systems. A combination of the present 
method with electronic transport and electron-phonon 
interactions is worthy of further explorations. 
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